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Extending the transfer matrix DMRG algorithm, we are able to calculate imaginary time spin 
■ autocorrelations with high accuracy (absolute error < 10~ 6 ) over a wide temperature range (0 < 

p. J < 20). After analytic continuation using the rules of probability theory along with the entropic 
prior (MaxEnt), we obtain real frequency spectra for the XY model, the isotropic Heisenberg and the 
gaped Heisenberg-Ising model. Available exact results in some limits allow for a critical evaluation 
of the quality of answers expected from this procedure. We find that high precision data are still 
insufficient for resolving specific lineshapes such as low frequency divergences. However, the method 
is appropriate for identifying low temperature gaps and peak positions. 
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I. INTRODUCTION 



Finite temperature dynamic correlations are the link between experiment and theory. In spite of their relevance, it 
seems that for strongly interacting electronic or magnetic systems, no reliable, direct theoretical methods exist for their 
evaluation. This absence is particularly felt recently in the field of quasi one-dimensional magnetic materials, where 
excellent samples, detailed neutron scattering and NMR experimentscra call for a better understanding of dynamic 
form factors. In this context, numerical simulation techniques can provide valuable information. However, these are 
, also subject to rather severe limitations. On the one hand, real time methods based on the exact diagonalization of the 
Hamiltonian matrix have been restricted to small systems. As a consequence, the extrapolation to the thermodynamic 
limit is often unreliable,.-especially in the low temperature regime.u Related methods such as the moment expansion 
or the recursion methodu can provide reliable short time correlations, however, the extrapolation to long times is left 
, uncertain. Until now, real time methods have been insufficient for discussing low frequency properties, for instance, 
y— ( ■ they cannot decide about diffusive or ballistic transport £10 

On the other hand, imaginary time methods such as quantum Monte Carlou (QMC) have to face the ill-conditioned 
OS . analytical continuation problem. In this context, it is reasonable to expect that the higher accuracy of the transfer 



matrix DMRG calculationsJj'E_l free of stochastic exrocs, might help. This method has been applied successfully to 
several interesting (quasi-)one dimensional systemsJiJTij In this work we describe in detail an extension of the transfer 
matrix DMRG method proposed in Ref. [H] to calculate imaginary time correlationsEJ We show that these can be 
very accurately determined, especially for local correlations. In addition, we present an extended critical discussion on 
the potential of this method to obtain reliable real frequency spectra. Of course, the transfer matrix DMRG method 
q | is limited to quasi-one dimensional systems. 

In order to investigate whether the combination of transfer matrix DMRG and analytical continuation methods can 
provide accurate dynamic correlations, we focus on the s = 1/2 antiferromagnctic Heisenberg Hamiltonian: 



H = J^2(SfSf +1 + SfSf +1 + ASf 5f +1 ), (1) 



i 



where Sf = \af , <jf are the Pauli spin operators with components a = x, y, z at site I. In the following we will take 
J as the unit of energy. 

This model represents a good playground for tests since in the XY limit (A = 0) exact results at all temperatures 
for the longitudinal zz-correlation and at T = 0, oo for the transverse xx-correlation are known. For A — 1, and A > 1 
where therSnectrum is gaped, the two spinon contribution to the transverse correlation at T — was recently exactly 
evaluated .EIE3 

The paper is organized as follows: section [n] explains how to extend the transfer matrix DMRG technique for 



obta ining imaginary time correlations and also provides a summary of the analytical continuation methods. In [II A 



and III B, we look closely at the XY model and compare our results to exact solutions. Among various analytical 
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continuation procedures, we fin d the MaxEnt method to be the most reliable. The resu lts for the isotropic point 
(A = 1) are presented in section HI C and those for the gaped regime (A = 2,4) in section III D . 



II. TECHNIQUE 
A. Notation 

Before explaining the transfer matrix DMRG technique, let us first fix some notation. The basic quantity from which 
the nuclear spin relaxation rate 1/Ti and the neutron scattering cross section can be determined is the dynamical 
structure factor 

/+oo 
e**(S?(t)S?(p))*t (2) 
-oo 



where a = z for the longitudinal and a = x for the transverse correlations. S"(t) — e ltH Sfe~ ttH and the average 
(•) is taken in the canonical ensemble at the inverse temperature (3 = 1/T. S°j(u) is a positive function and the 
autocorrelations i = j satisfy the sum-rule 
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at all temperatures. Using the transfer matrix DMRG method we will study the imaginary time Green's function 

G%{t) = (S«{t)S«) (4) 

where Sf (r) = e TH Sfe- TH . 

It is related to the real frequency correlations through the linear integral equation: 

G U T ) = ^Z r K(T,u)S%(u)du (5) 
Z7r Jo 

K(t,lu) = e -T« + e -(fl-r)w ( 6 ) 

where the detailed balance condition S^{— w) = e~^ u Sf^(u>) has been included in the kernel K(t,uj). Notice also 
that using the symmetry K{t,uS] — Kip — t,ui) it is only necessary to calculate G?At) in the interval [0,/3/2]. The 
decomposition of the integrand in Eq. into K(t, uj) and Sfj (oj) is not unique. In fact, it is sometimes better (see the 
xx correlations in the XY model) to reconstruct the symmetric correlation function Sfj{u), even in u), corresponding 
to the Fourier transform of the anti- commutator correlations. For this purpose, we rewrite Eq. (^|) as 
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G U T ) = ^ll K{T^)S^)du (7) 



o 



with 



k(T,u) = K(T,cj)/(l + e-P"), 

Sg.( W ) = SgH(l + e-^). (8) 



This scheme will be referred to as the symmetric scheme. 



B. Transfer matrix DMRG 



In this purely technical part, we will explain in detail how one can calculate imaginary time correlations by extending 
the transfer matrix DMRG method developed in reference [ [lC}]. 
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1 . Thermodynamics 



Let us briefly recall how the transfer matrix representation of the partition function is used to calculate thermody- 
namic quantities. First, we split the Hamiltonian (|l]) into odd and even bond terms: 

H = H Q + H e 

H a = h 1 + h 3 + h 5 + --- 

H e = h 2 + Ha + h 6 + ■ ■ ■ 

hi = S*S* +1 + SVSf +1 + AS*S* +1 , 

so that the all even (respectively odd) terms hi commute with each other. Then, using the Trotter-Suzuki decom- 
position and doing the standard insertion of the identity 1 = Y\i=i } l s fc')( s fc'l & t each inverse temperature slice 
k' = 1, 2, . . . , 2M, we can express the partition function in terms of the quantum transfer matrix Tm- 

Z = Tre-? 11 = Tr [e~ eH "e-^] M + 0(e 2 ) , 

Tr[e-^e-^] M = Tr [T M ] N/2 - (9) 

Here M is the Trotter number, e = (3/M and j3 = 1/T. In the last step of Eq.(||), the summation indices have been 
permuted so that the space and Trotter (inverse temperature) directions are interchanged. 7m is a non-symmetric 
matrix given by the product of 2M local transfer matrices 

M 

(°1 • • •°'2mI^m| '1 • • -°2Af) = ^2 II X ( CT 2fc-l CT 2fcl CT 2fc-l CT lfc)^ r ( cr 2fc cr 2fc+ll cr 2fc cr lfc+l) ( 10 ) 

with periodic boundary condition in the Trotter direction ((72M+1 = fi), cr^. = (— l) t+fc s' fc and 

r(44 + xWi +1 4 + +\) = (4+i»4 + i|exp(-e/ ii )|4,4 +1 ). (11) 

The real space indices are denoted by i, the Trotter direction by k and Sf \ si ) = s l k \s l k ). The sign in the definition of 
a l k is chosen such that r conserves er: a k +cTfe +1 = v k +1 +^k+i- Eq s - (|9|)-(|rT|) are best represented as a 2-dimensional 
checkerboard shown in Fig. |l|(a) . The arrow in the left corner emphasizes that r propagates in the real space direction. 
Fig. |](b) shows how the full checkerboard is cut along the real space direction in N/2 identical quantum transfer 
matrices Tm- 

In the limit N — > oo, the partition function Z is given by the maximum eigenvalue Amax- Thermodynamic quantities 
such as the magnetization or the internal energy can be obtained from the corresponding left (ip L \ and right \ip R ) 
eigenvectors of the transfer matrix Tj,/ .t3 

Let us emphasize that the method provides results free of finite size errors, the thermodynamic limit being obtained 
automatically, due to the exponent N/2 in Eq. (||). Remaining errors are the systematic 0(e 2 ) Trotter error and the 
systematic error introduced by the truncation of the basis set for the density matrix.c3 

2. Imaginary time correlation function 

We now turn to the imaginary time correlation function (|J) and restrict ourselves to (r) , the extension to Gfj (r) 
being straightforward. For convenience, we express (Q) as 

OUr) = ^ 

Mj(t) = Tr (s*e- TH S*e-^- T)H ^j = Tr (S*(e- eH ) k S*(e~ eH ) M - k ) , (12) 
where r = ek and fc = 0, 1, . . . , M — 1. Using the following symmetric decomposition: 

e~ cH = e~i H °e~ eH e e -i H ° +0(e 2 ), (13) 

we have 
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Afij (r) = Tr ( e-* g °Sfe-* g °e- 6g ° e- gg °e- gg ', • ■ ■ er tH °er tH \ 

l 2 k 

e - iHoS z e -jH ee - e H. - 6 g. . , , e-*H -eH. j _ (w) 

up to C(e 2 ) corrections. 

As for the partition function, we insert identities along the horizontal Trotter slices, permute the summation indices 
in the trace and form local transfer matrix chains along the Trotter direction to obtain an expression similar to Eq. (0) : 

A/-,(r) = TrK|(fc) (T M ) ( £-[iMlH)]. (15) 

Here [§] is the closest integer larger than or equal to |. In defining T^f (k), we have to differentiate the case j = i, i + 1 
from j > i + 1. First, however, we need to define the local spin transfer matrices 

n z (4,4 + M + \4 + +\) = <4>4 + V f ^f e -^|4 +1 ,4+\) (ie) 

^(4,4+il4 +1 .4+i) = (4, si +1 \e-i h *S? S*e-4*'|4+i>4+i> (17) 
where I, m = i, i + 1. Then, for j > i + 1, 

r^(fc) = r A v(o) (T M )[iHiH 7&(*o (is) 

with 

T H a 2fc+lJ CT 2fc+2l CT 2fe+l' CT 2fc+2J T V a 2fc+2' a 2fe+3l a 2/c+2' ^fe+S/ 1 

X 11 T (CT2fe'-1^2 fc '|CT 2 fc'-l' CT 2fcO T ( cr 2fc' '^fc' + ll^fc' '^fe' + l)' 
fe'G[l,Af] 

Again, these formulas can be nicely understood graphically as in Figs. 0(c) and 0(d), corresponding to Eqs. 

for the case N = 6, M = 3,i = l and j — 3,4. In this case, the two local spin transfer matrices sit on different 

transfer matrix chains. 

On the contrary, for j = i, i + 1, they belong to the same chain: 

K--^2Ml^(fc)kl +2 --^2M) (20) 



T i ("l.^K >°2 M°2 >°3 l°2 >°3 

{-; +1 i 

V t Z ('^r* /r' \ |_i+2 ^r* + 2 \ 

x T ml%+1' CT 2fc+2l cr 2fc+l' (J 2fe+2- ,r| > Cr 2fc+2' Cr 2fc+3l Cr 2fc+2' "2k+3> 
k'^l.k+l 

X 11 r ( Cr 2fe'-l> Cr 2fc'l Cr 2/c'-l' (J 2l' M^fc' > CT 2fc'+ll°2fc' > °2fc'+l) 



fe'G[l,M] 



when fe = 1, 2, . . . , M — 1 and 



K--^2Ml^(fc)K +2 --^2M) (21) 

= E ^L(4,4\4 + \4 +1 M4 + \4 +1 \4 +2 ,4 +2 ) 

X 11 1^2^-1' 2fe' I 2fe'-l' Sfe' ) ( 2fc' ' 2fe'+ll 2fc' ' 2fc' + l)' 

fc'e[i,M] 



for fc = 0, corresponding to a static correlation. 

Finally, we obtain for the imaginary time correlation 
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G?-(r) 

In the limits N — ► oo and |i — j\ <C AT, this reduces to 



Tr[T^(fc) (7 M )(f-[l] + W- 1 ) 



Tr(T, 



(^ir^(fc)i^) 
,m-m+i 



(22) 



(23) 



For systems with finite correlation length £, Gf.{r) ~ (S?) when |j — i| 3> £, which can be verified systematically 
as £ is determined from £ — 1 = i In 



A„ is the next largest eigenvalue of Tm- 



3. Renormalization of transfer matrices 



Now that we have defined the relevant transfer matrices, we must explain how to construct them in practise. The 
purpose is to add successively new r plaquettes to grow Tm in the Txotter direction (in the same spirit real sites 
are added to the Hamiltonian in the standard T = DMRG algorithm^) and simultaneously find a renormalization 
procedure that keeps the dimension of the matrix 7u fixed as the temperature is lowered. For this purpose, we cut 
Tm schematically in two halves as shown in Fig. || for the two generic cases of M odd (a) and even (b). In the 
DMRG language, Tm is called the superblock, the right inner block (dashed line in Fig.||) plus the left edge (<Ti, a 1 ) 
the system and the left inner block plus the right edge (02, <x 2 ) the environment. We use n s and n e to label the basis 
sets in the inner of the system and environment blocks, respectively. The states at the left and right edges are labelled 
by a\ and 02- With this notation, we denote the elements of the right transfer matrix by S (a'{, n' B , cr'2', <J\, n s , 02) or 
S e (a{,n' s , cr 2 '; a", n s ,02) depending on whether the system consist of an even (e) or odd (o) number of r plaquettes. 
When adding a new plaquette, the elements of the new system are given by the following recursion relation: 

5 e (ffi, h' s , 02 5 a",n s ,a 2 )=J2 T( t r' 1 ,o'\o'{,o")S (o",n'„ a, n s , <r 2 ), (24) 

a" 

S (a'{, n' s , cr'2] <ri,n a , a 2) = ^ T K. (T 'Vi>' T ) 5 e(' T '. n' s , of; a", n s ,a 2 ), (25) 

where {\n s )} = {\a}} <g> {\n s )}. Initially, when AS = 2, S e (a[, a' , a%; of , a, a 2 ) = T{a , l ,a , \a , {,a ,, )T{a ,, ,a^\a 1 a a ). 
Similarly, the elements of the left transfer matrix £ (a' 1 , n' e , a 2 ; u'{ , n e , a'2) or E e {a'{, n' e , a^; <J\, n e , cr 2 '), are given by 

£ e , K , <J 2 ; <Tl , fl e , a 2) = ^ T (°1 ! G ' ' I °1 > W ! n 'e i °2 5 °"" > "e ) °2 ) : ( 26 ) 

£ K . n' e , o- 2 ; a'l , n e , a' 2 ' ) = ^ r fa , I ^ , <?")£e {a" , n' e , a' 2 ; <t, n e , cx 2 ') (27) 

<T /y 

after one plaquette has been added, {|n e }} = <S> {|n e )}. 

When the number of states in {|n s }} ({|^ e }}) exceeds a given number m, the transfer matrices (f e .o) are 
renormalized by: 

•4e(oi, <, <J 2 ;cr'{,n s ,a 2 ) = 0^(<,n'JAW,n' s ,CT 2 ';cri',n s ,CT2)0^(n s ,n s ), (28) 
^Kii's^aVi,".!^) = ^2 0\{n' s ,n' s )A {cr'l,h' s ,a'2^cr 1 ,n s ,a2)O r A {n s ,n s ), (29) 

.4 = S for the system block or £ for the environment block, n' s , n s — 1, 2, . . . , m. The transformation matrices O l A and 
are constructed from the m largest eigenvectors of the corresponding reduced, non-symmetric density matrices 
given in terms of the left and right eigenvectors of Tm by: 

p s =Tr„ c ^ 2 |V fl ')(^ L |, Pe =Tr na , (71 |^ R )(^ i |, (30) 

for the system and the environment. In terms of S and £ , the renormalized transfer matrix Tm corresponding to the 
superblock (Fig. ||) is given by 
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£ {a l l ,n' e ,a l 2 ;a'l,n e ,a'.i) 



TM(n' e ,<T2,(r[,ri' s ;ne,<T2,(Ti,n s ) = < 



x S e (a[,n l s ,a 2 l ;a l {,n s ,a 2 ) 



(31) 



for M/2 odd or even. 

Many systems of interest have spatial reflection symmetry. For these, we can obtain the left transfer matrix from 
transposition of the right one so that for instance, £ {o\, n',a' 2 ; <t", n e , Ojf) = S (a'l, n e , a 2 \ ax,n' e , a 2 )- Consequently, 
the left eigenvector of (^ L | can be obtained from \ip R ): ip L (n s ,<j2,o'i,n>e) — ^ R (n e , &2, &i, n s ). For systems without 
reflection symmetry such as zigzag or dimerized chains, left and right transfer matrices and eigenvectors must be 
evaluated separately. 

As a last technical step, we have to explain how to renormalize T^j(k), needed for evaluating G?-(r). Again, we 
must distinguish two cases: 

For 3 = i, i + 1, the two local spin transfer matrices are located on the same transfer matrix chain. In this case, we 
can separate T^(k) into left £^ {k) and right S\ Q parts, similar to what we did for 7m- As G^(r) = G^ifi — r), it is 

sufficient if k = 0, 1, . . . , M/2 so that there is always exactly one local spin transfer matrix in each half of T^j (k) (Fig. 
§.) The blocks £{ (k) and S l e o are renormalized according to Eqs. (^5|)-(^9|) except that t has to be substituted with 
rf at the appropriate steps. The static, k = 0, case has to be treated separately. 

When j > i + 1, T^|(fc) consists of ([fWfl+l) adjacent parallel chains which must be renormalized as a whole. 
Again, they can be cut into halves £ 3 e {k) and S l e a , with ni = 2([§]-[|]) + 1 internal real space lines instead of one 
in the previous case (Fig. ^J). These can be renormalized according to Eqs. (|25|)-(|2Sj) with the proper modifications. 
However, the storage of £l a (k), S l e will be increased by a factor (25 + Xy( n ^*) which will eventually become 
restrictive for large j — i. If this is the case, we approximate G^(r) by substituting the renormalized Tu, 7^(0) and 

T^(k) in the expression for T^(k) (Eq. jl^)) and multiply them successively to \ip R ) in Eq. (|23|). The accuracy of 
the obtained result can be controlled by varying the number of states m. 

For a realistic description of NMR experiments, which involves at most j = i, i+ 1, i + 2 correlations, the calculation 
can be easily performed for systems having up to three states per site. 



C. Analytical continuation 

The analytical continuation is concerned with the inversion of the integral equation (Q) whose principal feature is 
that the kernel K(t,u>) is very singular. This renders a numerical inversion particularly sensible to errors in Gfj(r) 
which are exponentially amplified in the result. Before discussing the probabilistic approach to ill-posed inversion 
problems we will briefly discuss the SVD method because it gives insight on how ill-posed the problem is. 

Let us first discretize the r and lo variable and restate Eq. (|^) in matrix form, omitting the superfluous subscripts 
for now, 

G(r i )=^2K ij S(u; j ), i = l,...,N, j = l,...,M (32) 

j 

At this stage we should mention that it is straightforward to incorporate knowledge on the derivatives of G^ (r) = 
■jptG(t) by just adding more equations 

G ( "V0 = E A lT )5 K)' l=l,-,N', j = l,... ,M (33) 

j 

to the linear system (^2|). The formal problem is unchanged, only the vector of data points and the kernel are larger 
in the first index. 

The SVD decomposition of the matrix K is = 53/=i UuAtVji where V is an M x M orthogonal matrix, while 
the N x M matrix U is merely column-orthonormal, i.e. U T U is the M x M unit matrix. The Af correspond to the 
M eigenvalues of the matrix K T K. Formally, the solution of equation (|32|) is: 

M N 

s(^) = Y,J2 v ^Y UaG{n) (34) 

1 = 1 i=l 1 



G 



One immediately anticipates the catastrophe when the eigenvalues A; become small and G(r) contains errors. Since 
Kij > 0, small eigenvalues correspond to rapidly oscillating eigenvectors of K T K which therefore couple strongly to 
the noise. To illustrate the ill-posed nature of the analytic continuations, consider the situation where N = M = 100. 
The condition number , i.e. the ratio of the largest to the smallest |Aj|, is greater than 10 17 for the entire range 
of interesting /3-values . In other words, even the errors introduced by the finite machine-accuracy are sufficient to 
make the direct inversion useless. For f3 = 1 only 5% of the eigenvalues are greater than 10 -8 and for (3 — 10 the 
percentage is 10%. The situation is roughly unchanged when including derivatives. The straightforward application 
of Eq. ([m]) yields results which are orders of magnitude too large. The natural way to regularize the sum (Eq. |34| ) 
is by truncating it at l cut so that the error 5(J2i G( T i)Uu eut ) <C Yli G{Ti)Uu cut . Typical values of l cut for the transfer 
matrix DMRG data are 7 to 10. The drawback of this ad-hoc truncation-scheme is that a major part of the vector 
space for S(ojj) is lost. A further disadvantage of the SVD-approach is that it does not enforce positivity of the 
solutions. It should be noted that the SVD-approach is equivalent to Tichonov-regularization, where the L 2 -norm of 
the image is the regularization criterion. We should also mention that we tested the Pade approximation. It appeared 
that it did well in some situations (isotropic case), but it failed in others (gaped phase). Therefore, we will not discuss 
this method further. 

There is a wide-spread misconception about ill-posed inversion problems, or inductive inference problems in general. 
It makes no sense to ask for the true function f(u). There is no chance, whatsoever, to infer the true result! All 
inference schemes yield results which can in principle deviate widely from the unknown true result. The correct 
question to be asked and which can uniquely be answered is rather: what is the distribution of functions f{to) 
compatible with the noisy and incomplete data and all our prior knowledge. In other words, we should aim for the 
probability density p(f(uj)\D, I) for a function f{u>) in the light of the transfer matrix DMRG data D and additional 
prior-knowledge /, such as sum-rules and positivity constraints. The elementary product-rule of probability theory 
allows us to determine this probability consistently: 

p(f(u 1 )\D,I)=p(f(u 1 )\I)p(D\f(u),I)/p(D\I) (35) 

in terms of the likelihood p(D\f(ui), I), the prior p(f(u>)\I) and the normalization p(D\I). The likelihood stands for the 
probability for the data D, assuming that f{ui) is the exact function. The likelihood deviates from a delta- functional if 
the data suffer from statistical noise or unknown systematic errors, like in the present case. For the likelihood the source 
of the missing information does not matter. As long as we know nothing about the features of the systematic errnts, 
the data have to be considered as the mean, and the errors as the variance of the likelihood-distribution functioned. 
Hence, like in the case of uncorrelated normal-distributed data, the likelihood reads p(D\f(ui), I) cx exp(— x 2 /2), with 

2 ^ |G(n)-E 3 ^MI 2 

x =2_ -2 ■ ( 36 ) 

1 * 

The prior distribution p(f(to)\I) quantifies our knowledge about the solution / prior to our knowing the data D. For 



a general scheme, the only reliable prior knowledge is that / is a positive 
adequate probability distribution is the entropic prior p{f{u>)) oc exp(aS). 
or relative entropy 



27 



ditive distribution function for which the 
with S being the information divergence 
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/U',1- mu,) (^gy) ) A',,, C!7) 



of the function f{oJj) relative to a default model m(ujj). The factor Koj accounts for the negative frequency contribu- 
tion consistently with the detailed balance condition. Maximizing the posterior probability p(f {lj)\D , /) is equivalent 
to maximizing the functional 

cj> = aS - l - X 2 (38) 

When maximizing </>, S essentially regularizes the solution such that: (i) it stays positive and, (ii) structure relative to 
the default model is penalized, depending on the parameter a. MaxEnt yields the most uncommittal solution which 
shows only structure if it is significantly supported by the data. In "classic" MaxEnt, a is determined self-consistently 
using the rules of probability theory, i.e. such that the solution f(u)j) is the most probable in the light of the input 
data (details can be found in Ref . [ |29| , |30[ | ) . We have only considered a flat m{u>j) for loj > 0. 

As mentioned above, in the transfer matrix calculation (as opposed to QMC), the errors cr, are systematic but 
unknown. They are due to the truncation of the basis set and to a finite Trotter step. Since the values for oi are 
not known we have to determine the probability p((Ji\D,I) for Oi using the rules of probability theorya. Typically, 
<7i ~ 10 -6 for (3 — 16. So doing, the reconstructed imaginary time correlation agrees with the DMRG data up to 
- 10~ 6 or better. 
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III. RESULTS 



A. XY model: Longitudinal autocorrelation 

The XY model (A = 0), which can be mapped onto a free fermion model via a Jordan- Wigaer transformation, is 
useful for tests because its longitudinal zz— correlations in (q, uj) can be expressed in closed formE3 at any temperature 
T. 

The corresponding imaginary time autocorrelation function can be represented as 

(1 rcosq \ 2 

- / z dq) (39) 

In Fig. ||, we compare this exact result with the transfer matrix DMRG data for (3 — 2, 8, 20 on a logarithmic 
scale. The reversed peaks are merely artifacts due to the change of sign in the argument of the logarithm. For all the 
calculations, we have kept m — 100 states in the density matrix so that the truncation error (1 — Y] m p m ) is smaller 
than 10 -7 for the largest Trotter number M=800 ((3 = 20 when e = 0.025). In order to reduce the Trotter error, we 
have done a linear e 2 — > extrapolation (this requires commensurate values of e). This procedure is justified as long 
as the systematic errors induced by the truncation are negligible compared to the Trotter errors. This is typically the 
case at high temperatures (small truncation error, Fig. ||(d)) or when e is large (large Trotter error, Fig. ||(a)). Fig. 
||(b) is an example where the e 2 extrapolation fails because the systematic errors become comparable to the Trotter 
error itself (M = 800 Trotter steps are needed to reach (3 — 20 when e = 0.025). In such a case, nothing can be gained 
from the extrapolation and it is justified to use the smallest e data for the analytical continuation. Notice that the 
result in Fig. g(a) obtained from the e = 0.1, 0.2 data is as precise as the e = 0.025 calculation. This can be exploited 
for models with more degrees of freedom per site, where small e calculations cannot be afforded and values of e = 0.1 
or larger may be needed to reach the low temperature regime. 

Before doing the continuation to real frequencies, we should warn that because the XY model essentially describes 
free fermions, its zz— correlation function is not generic of a true interacting system and exhibits some peculiar 
behavior. For instance, the zz— autocorrelation has a sharp steplike cutoff at uj = 2 and a uj — logarithmic 
divergence for (3 = 0. As we are not interested in reproducing step discontinuities, not expected in interacting systems 
at finite temperature, we artificially cutoff the integral at 10 — 2 in Eq. (|^) to investigate how the smooth lineshape 
in the interval uj £ [0, 2] can be reconstructed. 

In Fig. @, we compare the MaxEnt results obtained from the transfer matrix DMRG data, those obtained from the 
exact G(t) and the exact Sfl(uj). Notice that MaxEnt cannot resolve the low frequency divergence at (3 = 2 since 
there is not enough spectral weight under the logarithmic divergence. The MaxEnt results from the numerical data 
and those from the exact G(r) cannot be distinguished for (3 — 2 and (3 = 8. 

B. XY model: Transverse autocorrelation 

After the Jordan- Wagner transformation, the xir-correlations become nonlocal fermion correlations which are ex- 
pressible as PfaffiansEy at all temperatures. Closed form expressions are known only at T = 0, ooH The exact T = 
autocorrelations have singular behavior at all integer frequencies uj = I, I = 0,1,2,.... Sf^uj) diverges as uj^ 1 ! 2 for 
uj — > + and as In \uj — 1 for uj — > 1. The higher singularities are cusps, for uj — » 2 + , Sf^uj) ~ (uj — 2) 1 ' 2 . At T = oo, 

In Fig. uy we show the numerical data continued with MaxEnt in the intermediate to low temperature range 
(3 = 6 — 16. As the temperature is lowered, one can see the emergence of two peaks, one at uj = 1 and the other near 
uj = 0. Also, we see some sign of the cusp at u> = 2. Reducing the temperature, the peak near uj = continuously 
grows and shifts closer to the zero temperature singularity. That at uj = 1 consistently moves towards the correct 
T = position, as long as the temperature is not too low ((3 < 16). When the quality of the data becomes poor, 
which can be due to a large Trotter number M (low temperature), to a small basis set (small m), or to a large Trotter 
step e, we observe a systematic shift of high frequency structures towards lower frequencies. This effect can be seen 
in the inset when only m = 100 states are kept instead of m = 160 and e = 0.05 rather than 0.025. |— . 

Further, we learned that it is better to reconstruct Sf^uj) using the symmetric scheme. Indeed, it is knowna that 
the time correlation Sfj(t) decays exponentially for all and temperatures T, therefore, Sf^uj) is regular at uj = 0. 

This implies that ^5*^(0) = ^Sf { (0) due to the detailed balance condition. When reconstructing directly Sf^uj), this 
relation is not fulfilled (the two procedures are compared in the inset for [3 = 16). 
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Finally, in an attempt to increase the information contained in the input data, we calculated -^G{t) with equal 
precision as G(t). We find that the additional data does not help much in the analytical continuation (which can 
be understood when looking at the SVD decomposition of the extended kernel), except that it gives a hint that the 
curve obtained by the symmetric kernel is better behaved near us = 0. In any case, we found it more economical to 
improve the accuracy on G(r) rather than calculate -xpG(r). 

C. Isotropic Heisenberg model 

Rigorous results for the dynamic correlations in the isotropic Heisenberg model (A = 1) are rare. Recently^] the 
exact two spinon contribution to S x (q, us) at T = was found. In this restricted subspace, the autocorrelation function 
has singularities ln(w) for us — > and (In \ir/2 — us\) 3 ^ 2 as us ^ n/2. 

In Fig. |[ we present results for SJ-(w)- The main peak at us = tt/2 for T = is shifted to slightly lower frequencies 
at finite temperatures. As in the xx— correlations of the XY— model, this peak seems to move away from the exact 
value for the lowest temperatures studied, presumably due to loss of accuracy. Furthermore, in accord with the 
observation by Starykh et alia using the QMC method (/? = 8), a low frequency secondary peak develops when the 
temperature is lowered. 

We should also observe that the way the low frequency peak develops is much different from the zee-correlation 
in the XY case. There, the low frequency maximum can be traced to the fact that we chose the non-symmetric 
representation (Eq. (^)). Indeed, as the symmetric spectra have a maximum at us = for all temperatures, division 
by (1 + e^' 3 ") results in a low frequency maximum in the non-symmetric correlation. In contrast, in the A = 1 
model, the symmetric Sf^us) has a minimum at us — and a small low-frequency peak which survives upon division 
by (l + e"' 3 "). 

The zero frequency limit Sf^us — > 0) = Sf^us — > 0), relevant for NMR experiments, seems to be roughly temperature 
independent. In contrast, it increases sensibly with decreasing temperature in the XY— model. 

D. Heisenberg-Ising model 

The Heisenberg-Ising model (A > 1) iSrGharacterized by a gap in the excitation spectrum with values E gap = 0.39 J 
for A = 2 and E gap = 2.15J for A = 4.E3 At T = 0, both Sf^us) and S*(u) will exhibit the gap, however, a 6(us) 
function will subsist in S^(uj) due to the non- vanishing matrix element between the two degenerate ground states 
in the thermodynamic limit with total momentum quantum number k = 0, it. This matrix element was evaluated 

exactly by BaxterS (k = 0\Sf\k = n) = §n£U (tT^") where A = 2± §— • We have verified that G z u (/3/2) 
approaches \{k — 0\S?\k — ir)\ 2 at low temperature. When reconstructing S^us), the large weight 5(us) function 
renders the resolution at finite frequencies poor so that the gap cannot be seen as sharply as in Sf^us) shown in Fig. 

. Although the value of the gap can be estimated from the spectra, the precision is not comparable to the zero 
temperature DMRG method. We should also mention that thej-spectrum above the gap edge seems to start without 
discontinuity, as is suggested from the two spinon contribution.!^ This behavior helps in the identification of the gap 
value, in contrast to other models which exhibit singularities at the gap edge at T — 0. For instance, the S = 1 
chain shows a square root divergence, when assuming a constant matrix element around the quadratic minimum of 
the lowest magnon excitations at q — tt. Such divergences are smeared out in the MaxEnt analysis. 

IV. CONCLUSION 

We have investigated in detail whether high accuracy imaginary time data, obtained from the transfer matrix 
DMRG method, can be exploited in evaluating real frequency correlations after an analytical continuation. Using 
as a test model the the spin 1/2 Heisenberg chain and the MaxEnt method, we found using the SVD and MaxEnt 
methods, we found that features such as the location of peaks and gaps can be reliably determined. More quantitative 
information about precise lineshapes or the nature of divergences seems to lay beyond this procedure. 

Reliable real time methods need to be developed in order to overcome the severe intrinsic limitations of imaginary 
time calculations. 
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FIG. 1. The checkerboard representation for M = 3, N = 6. (a) and (b) represent Eqs. ([]) and (0), (c) and (d) Eqs. 
and © for i = 1, j = 3, 4. 
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FIG. 2. The superblock Tm is cut in a system S and environment £ block for M=5 odd (a), M=6 even (b). 
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FIG. 8. Isotropic Heisenberg Model. MaxEnt from numerical data with m = 80 and e — 0.025. 
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FIG. 9. Anisotropic model, A = 2, 4. MaxEnt from numerical data with m = 100 and e = 0.05. 
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